#! /usr/bin/perl -w

######################################################################
# Time-stamp: <04/08/19 13:20:09 ostolop>
# used for Package: EP::EPC::Core::Threeway_Similarity_Analysis
# Sequence Alignment Similarity Precomputation 
######################################################################


use strict;

package Bio::DB::GFF::Adaptor::faux;

use Bio::DB::GFF::Adaptor::dbi::pg;
use vars '@ISA';
@ISA = 'Bio::DB::GFF:Adaptor::dbi::pg';


package main;

use lib '/usr/local/EPNG/cgi-lib/perl/';
use EP::Config;

use DBI;
use Bio::Factory::EMBOSS;
use Bio::Tools::Run::EMBOSSApplication;

use SVG;

# constant
use constant MY_DEBUG => 0;

# constant for GFF database
use constant PG_DBI_DRIVER => "dbi:Pg:dbname=";
use constant PG_ADAPTOR => "dbi::Pg";
use constant DBNAME => "gff";
use constant USER => "postgres";
use constant PWD => "";

# constant for alignment
use constant ALIGN_FILE_EXT => ".aln";
use constant SUMMARY_FILE_EXT => ".summary";
use constant BASE_ALIGN_DIR => $EP::Config::h->{TSA_DATA_PATH} . "/align/";


#------------------------------------------------------------------
# 1. Compute the sequence alignment for all pairs of genes (GFF DB)  
#    and generate the pairwise alignment file for them
#------------------------------------------------------------------
    &all_to_all_pairwise_alignment();

#---------------------------------------------------------
# 2. Generate the summary similarity file for each gene  
#---------------------------------------------------------
    &batch_align_info_summary(BASE_ALIGN_DIR); 

#---------------------------------------------------------
# 3. GENERATE SEQUENCE ALIGNMENT SVG FILES 
#---------------------------------------------------------
    &precompute_seq_align_svg(BASE_ALIGN_DIR);


#################################################################################
# TimeStamp 2004-08-03 Phil Hyoun Lee
#
# METHOD     batch_align_info_summary      
# PURPOSE    make the ordered align score summary file for all gene 
# ARGS	     $dir_path
#            ex) $dir_path="/home/my_perl/IGA/align/";
# RETURN     
##################################################################################
sub batch_align_info_summary {

        print "\n\n********************************************\n";
        print "batch_align_info_summary() ....\n";
        print "********************************************\n";

	my $dir_path = shift;

	my ($gname_path, $output_path, $filename, $valid);
	
	# get the directory information from $dir_path
	opendir DH, $dir_path 
		or die "[ERROR] Cannot open the dir $dir_path : $!";

	foreach $filename (readdir DH) {
		# check if it is the directory
		next if -f $filename; 

		# check the directory is the gname path
		if( $filename =~ /^Y/ ) {    # could be changed 
			$gname_path = $dir_path . $filename . "/";
			$output_path = $dir_path . $filename . SUMMARY_FILE_EXT;

			&make_align_info_summary_for_one_gene($gname_path, $output_path);
 		
			print "$gname_path\n$output_path\n"; 
		}
	}

	closedir DH;
}     



#################################################################################
# TimeStamp 2004-08-03 Phil Hyoun Lee
#
# METHOD     align_info_summary_for_one_gene      
# PURPOSE    make the ordered align score summary file for one gene 
# ARGS	     $gname_path, $output_path, $align_file_ext   
#            ex) $gname_path="/home/my_perl/IGA/align/YAL069W/";
#                $output_path="/home/my_perl/IGA/align/YAL069W/YAL069W.summary";
# RETURN     
##################################################################################
my %summary;
sub init_summary {
	foreach my $key (keys %summary) {
		$summary{$key} = 0;
        }
}

sub make_align_info_summary_for_one_gene {

	my ($gname_path, $output_path) = (shift, shift);
        my @sorted_summary;

	my ($filename, $valid, %align_info, $pair_gname, $score, $max_score, $sorted_gname);

	# get the directory information from $gname_path
	opendir DH, $gname_path 
		or die "[ERROR] Cannot open the dir $gname_path : $!";

	# init the summary hash
	&init_summary();

	foreach $filename (readdir DH) {
		# check if it is the alignment file
		$valid = index($filename, ALIGN_FILE_EXT);
		if($valid>0) {
			%align_info = &parse_align_info_from_emboss_output_file($gname_path . $filename); 		

			# currently we use only score
			$pair_gname = $align_info{'2'};
			$score = $align_info{'Score'};
			
			# store at the summary hash object
			$summary{$pair_gname} = $score;

			if(MY_DEBUG) {
				print "$gname_path$filename\n";
				print "\t$pair_gname has score $score\n";
			}
			
		}
	}	

	closedir DH;


	# sort the summary according to the score decreasing order
	@sorted_summary = sort by_score keys %summary;

	if(MY_DEBUG) {
		foreach $sorted_gname (@sorted_summary) {
			print "\t$sorted_gname has score $summary{$sorted_gname}\n";
		}
	}

	# write the sorted summary on the output file
	open FH, ">$output_path" or
		die "[ERROR] cannot open the file " . $output_path . " : $!";

	$max_score = $summary{$sorted_summary[$#sorted_summary]};
	for( my $i=$#sorted_summary; $i>0; $i-- ) {
		$sorted_gname = $sorted_summary[$i];
		$score = $summary{$sorted_gname} / $max_score;

		# make the score as percentage.
		print FH "$sorted_gname\t$score\n";
	}

	close FH;

}     

sub by_score {
	$summary{$a} <=> $summary{$b}  
}





####################################################################
# TimeStamp 2004-08-03 Phil Hyoun Lee
#
# METHOD     parse_align_info_from_emboss_output_file      
# PURPOSE    parse EMBOSS pairwise alignment result from file  
# ARGS	     $filepath   ex. "/home/my_perl/IGA/align/YAL069W/YAL020C.aln"; 
# RETURN     %align_info : HASH   etc. Currently we use only Score 
####################################################################
sub parse_align_info_from_emboss_output_file {

	my $filepath = shift;
	my %align_info;
 
	#my FH;
	my ($line, $delimiter, $valid, $name, $value) = (undef, ":");
	

	# open the file 
        open FH, "<$filepath"  or 
		die "[ERROR] cannot open the file " . $filepath . " : $!";
	


	# read the content
	while( defined($line=<FH>) ) {
		
		chomp($line);

		# if this is the part which has the alignment info, parse it.
		if( $line =~ /^#/ ) {
			# check it has the parsing info. currently parse everything with colon.
			$valid = index($line, $delimiter);
			if($valid>0) {
				# parsing start!
				# remove the comment notation and separate it
				$name = substr( $line, 2, $valid-2 );	
				$value = substr( $line, $valid+2, length($line)-$valid );

				# store it at return value
				$align_info{$name} = $value;
				
				if(MY_DEBUG) { 
					print $name, "=", $value, "\n";
				}
			}
		}
		# when it meet the first alingment sequence, stop this process.
		elsif( length(chomp($line))>5 ) {   # 5 is determined arbitrary !
		    last;
		}
		
	}	

	# close the file
	close FH;	

	if(MY_DEBUG) { 
		while( ($name, $value) = each %align_info ) {
			print $name, "=", $value, "\n";
		}
	}

	return %align_info;
}     


####################################################################
# TimeStamp 2004-08-02 Phil Hyoun Lee
#
# METHOD     all_to_all_pairwise_alignment      
# PURPOSE    batch pairwise alignment using EMBOSS package 
# ARGS	     $seqobj1 : Bio::Seq object  etc. should have -id, -seq 
# 	     $seqobj2 : Bio::Seq object  etc. should have -id, -seq 
# RETURN     $gid   
####################################################################
sub all_to_all_pairwise_alignment {

        print "\n\n********************************************\n";
        print "all_to_all_pairwise_alignment() ....\n";
        print "********************************************\n";

	my ($dbh, $sql, $sth, $gid, $row, $aligned_no, $result_no); 

	$dbh = DBI->connect( PG_DBI_DRIVER.DBNAME, 
                             USER, 
                             PWD,
		             { RaiseError=>1, AutoCommit=>0 } );

        #$sql = "select gid from fgroup where gname like \'Y\%\' and gid>=47 and gid<55";
        $sql = "select gid from fgroup where gname like \'Y\%\'";
	$sth = $dbh->prepare($sql);
	$sth->execute();
	
	$result_no=0;
	while($row = $sth->fetchrow_arrayref) {
		$gid = ${$row}[0];

		# show progress	
		print "GID " . $gid ;

	        $aligned_no = one_gene_to_all_pairwise_alignment_by_gid( $gid );
                if( $aligned_no > 0 ) {
		        print "\t success with aligned no " . $aligned_no . "\n";
			$result_no++; 
		}
		else {
		        print "\t fail with aligned no " . $aligned_no . "\n";
		}
	}

	# show result
	print "\n Total " . $result_no . " all_to_all processed \n";
	
	$dbh->disconnect;
}


####################################################################
# TimeStamp 2004-08-02 Phil Hyoun Lee
#
# METHOD     one_gene_to_all_pairwise_alignment_by_gid      
# PURPOSE    pairwise alignment using EMBOSS package 
# ARGS	     $gid :  
# RETURN        
####################################################################
sub one_gene_to_all_pairwise_alignment_by_gid {

	my $gid = shift;  
        my ($db, $gid_obj, $all_stream, $pair_obj, $pair_seq, $aligned_pair_no, $cmd);

	# connect to GFF DB
 	$db = Bio::DB::GFF->new( -adaptor=>PG_ADAPTOR, 
				 -dsn=>DBNAME, 
				 -user=>USER , 
				 -pass=>PWD )
	                    or die "Can't open database: " , Bio::DB::GFF->error, "\n";  

	# get sequence information from GFF DB 
	$gid_obj = $db->get_feature_by_gid( $gid ); 
	if(!defined($gid_obj)) {
		warn "[ERROR] gid ". $gid . " has no feature stream \n";
        }
	else {
		# everything is fine. make directory for storing output files
		$cmd = "mkdir -p " . BASE_ALIGN_DIR . $gid_obj->display_id ;
		system $cmd;

		$cmd = "chmod 777 " . BASE_ALIGN_DIR . $gid_obj->display_id ;
		system $cmd;
	}


	# get all Bio::seqIO streams from GFF DB
	$all_stream = $db->get_seq_stream( 'gene' );
	$aligned_pair_no=0;
	

	print "\nBase : $gid_obj computing for all other genes ..\n";

	while(1) {

		$pair_obj = $all_stream->next_seq;
		if(!defined($pair_obj)) {
			last;
		}  

                $pair_seq = $pair_obj->seq;


		# if all objects are processed, go out !


		# if it's not the proper gene, skip it !
		if(!( $pair_obj->display_id =~ /^Y/ )) { 
			next;
		}

		# if there is no sequence, something is wrong !
		if( length($pair_seq)==0 ) {
		        warn "[ERROR] gene " . $pair_obj->display_id . " has no sequence \n";
		}
		

		# everything is fine.
		if(MY_DEBUG) {
			print "Pair " . $pair_obj->display_id . 
                              " has sequence whose length is " . length($pair_seq) . "\n";
		}
		
		# this should work !
		#if( emboss_pairwise_alignment($gid_obj, $pair_obj) ) {
		emboss_pairwise_alignment($gid_obj, $pair_obj);
			# print $pair_obj->display_id , "\t";  
                	$aligned_pair_no++; 	

			# for drawing, store the sequence as a file.
			#my $seq_file_name = $id . ".seq";
			#open SEQFILE, ">".$seq_file_name or die "$!"; 
			#print SEQFILE $seq;			
			#close SEQFILE;

			#draw_alignment($seq_file_name);
	        #}

	}

	if(MY_DEBUG) {
		print "\n Total " . $aligned_pair_no . " pairs are aligned \n";
	}
	
	
	#$db->disconnect;

	return $aligned_pair_no;
}

####################################################################
# TimeStamp 2004-08-02 Phil Hyoun Lee
#
# METHOD     one_gene_to_all_pairwise_alignment_by_gname      
# PURPOSE    pairwise alignment using EMBOSS package 
# ARGS	     $gid :  
# RETURN        
####################################################################
sub one_gene_to_all_pairwise_alignment_by_gname {

	# input
	my $gname = @_;
	
        # get gene id from gene name
	my $gid = get_gid_from_gname( $gname );

        one_gene_to_all_pairwise_alignment_by_gid( $gid );

}


####################################################################
# TimeStamp 2004-08-02 Phil Hyoun Lee
# 
# METHOD     emboss_pairwise_alignment 
# PURPOSE    pairwise alignment using EMBOSS package 
# ARGS	     $seqobj1 : Bio::Seq object  etc. should have -id, -seq 
# 	     $seqobj2 : Bio::Seq object  etc. should have -id, -seq 
# RETURN     $result  : if program run successfully or not   
####################################################################
sub emboss_pairwise_alignment {

	# input
        my ($seqobj1, $seqobj2) = @_;

	# EMBOSS 
	my ($factory, $program, $program_outfile, $result); 

  	# get an EMBOSS application  object from the factory
  	$factory = Bio::Factory::EMBOSS->new(); 
  	$program = $factory->program('water');
        $program_outfile = BASE_ALIGN_DIR .  
                           $seqobj1->display_id."/".$seqobj2->display_id . ".aln";
  	if(MY_DEBUG) {
  		print $factory->program_info('water'), " program loaded !\n"; 
		print "input 1 : ", $seqobj1->display_id , "\n";
                print "input 2 : ", $seqobj2->display_id, "\n";
                print "output  : ", $program_outfile, "\n";
  	}

  	$result = $program->run({ -asequence => $seqobj1,
              -bsequence=> $seqobj2, 
	      -aformat => "pair",
 	      -awidth3 => 80,
              -gapopen   => '10.0',
              -gapextend => '0.5',
              -outfile   => $program_outfile});

	# now you might want to get the alignment, but useless !
        if(MY_DEBUG) {
  		use Bio::AlignIO;
  		my $alnin = new Bio::AlignIO( -format => 'emboss',
                                      -file   => $program_outfile );

  		while( my $aln = $alnin->next_aln ) {
               		# $aln is a Bio::SimpleAlign object
	       
               		print $aln->length, "\n";
               		print $aln->percentage_identity, "\n"; 
               		# this value is different with identity in output file !
               
	       		#print $aln->score, "\n";            # not initialized !
	       		#print $aln->no_residues, "\n";
               		#print $aln->is_flush, "\n";
               		#print $aln->no_sequences, "\n";
               		#print $aln->consensus_string(50), "\n";
  		}
	}
                        
	return ($result);
}

##########################################################
# TimeStamp 2004-08-02 Phil Hyoun Lee
# 
# PURPOSE    Get gene id from GFF database by gene name 
# ARGS	     $gname : ORF name	ex) 'YAL006C' 
# RETURN     $gid   
##########################################################
sub get_gid_from_gname {

	# input	
        my $gname = shift ;    
	my $gid;

	# DB 
	my ($dbh, $sql, $sth, $result_no, $row); 


	$dbh = DBI->connect( PG_DBI_DRIVER.DBNAME, 
                             USER, 
                             PWD,
		             { RaiseError=>1, AutoCommit=>0 } );

        $sql = "select gid from fgroup where gname=\'". $gname . "\'";
	$sth = $dbh->prepare($sql);
	$sth->execute();
	
	$result_no=0;
	while($row = $sth->fetchrow_arrayref) {
		$gid = ${$row}[0];
        	$result_no++;
	}

	$dbh->disconnect;

	if($result_no>1) {
		warn "[ERROR] Multiple gids retrieved from gname ",$gname, "\n";
	}


	if(MY_DEBUG) {
		print "get_gid_from_gname ", $gname, "=>", $gid, "\n\n"
	}

	return ($gid); 

}



#################################################################################
# TimeStamp 2004-08-09 Phil Hyoun Lee
#
# METHOD     precompute_seq_align_svg 
# PURPOSE    generate the sequence alignment svg files 
# ARGS	     $dir_path 
##################################################################################
sub precompute_seq_align_svg {

  print "\n\n********************************************\n";
  print "precompute_seq_align_svg() ....\n";
  print "********************************************\n";

  my $dir_path = shift;

  # get the directory information from $dir_path
  opendir DH, $dir_path 
		or die "[ERROR] Cannot open the dir $dir_path : $!";

  foreach my $filename (readdir DH) {

     # check if it is the directory
     next if $filename =~ /.summary$/; 

     # debugging 
     print "Directory : $dir_path" . "$filename\n";

     # check the directory is the gname path
     if( $filename =~ /^Y/ ) {    # could be changed 
	 my $sub_dir_path = $dir_path . $filename . "/";

         # make the temporary directory
         # system mkdir "$dir_path" . "$filename/";

         opendir SH, $sub_dir_path;
         foreach my $sub_filename (readdir SH) {
              # check if it is the directory
              if ( $sub_filename =~ /.aln$/ ) {   
                   # print "$sub_filename\n";
                   &generate_svg( $sub_dir_path . $sub_filename, 
                                  $filename . "/$sub_filename" );
              }
         }

         closedir SH;
     }
  }

  closedir DH;

}


#################################################################################
# TimeStamp 2004-08-09 Phil Hyoun Lee
#
# METHOD     generate_svg 
# PURPOSE    generate one sequence alignment svg file 
# ARGS	     $filepath : sequence alignment file
#            $out_filepath : svg output file
##################################################################################
sub generate_svg {

  my $filepath = shift;
  my $out_filepath = shift;

  my $svg = SVG->new( width => 500, height => 140, 'xml:space' => 'preserve', standalone => "no");
  my $line_no=0;
  my $text_height = 10;

  open FH, "<$filepath" or die "can't open the $filepath : $!";
  while(<FH>) {
      $line_no++;

      next if $line_no<25;
      last if $line_no>40;

      chomp;

      $svg->text( 
                  x => 0,
                  y => $text_height*($line_no-24),
                  style => "font-family: courier;" )->cdata("$_");
  }
  close FH;

  open FH, ">/usr/local/EPNG_DATA/TSA_DATA/align/$out_filepath.svg"  or die "can't open the $out_filepath : $!";
  print FH $svg->xmlify();
  close FH;

}
